/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
/** @file Arrays.cpp
 *  @brief Implementation of Vector and Matrix classes of double numbers
 */

#include "Arrays.h"
#include "Crout.h"
#include "Cholesky.h"

namespace smv_lib
{

//Defining ARRAYS_DEBUG turns on bounds checking and size checking
//#define ARRAYS_DEBUG
//================================================================= VECTOR CLASS
/**
** \brief class VECTOR of double numbers
**
**  VECTOR: VECTOR class of double numbers
**  Constructors: 
**  (1) empty (default) constructor: (size is assigned as default value 2, 
**      name is given as "VECTOR", but data should be given)
**  (2) constructor with a size 'n' and all elements 'd' (name is not given)
**  (3) constructor with a size 'n' and a given Name 
**  (4) constructor form a vector file (5) copy constructor
**/
//-------------------------------------------------------------------------
//
//! empty (default) constructor
Vector::Vector()
{
  data = new double[ M_DEFAULT ];
  assert( data != NULL );
  len = M_DEFAULT;
  Name = "VECTOR";
  isNameGiven = false;
}
//! Vector with 'n' double elements of value 'd'
Vector::Vector( int n, double d ) : len( n )
{
  data = new double[ n ];
  assert( data != NULL );
  for ( int i = 0; i < n; i++ )
    data[ i ] = d;
  Name = "VECTOR";
  isNameGiven = false;
}
//! Vector with given length 'n' and a name (all elements are set to zero)
Vector::Vector( int n, string name ) : len( n ), Name( name )
{
  data = new double[ n ];
  assert( data != NULL );
  double zero = 0.0;
  for ( int i = 0; i < n; i++ )
    data[ i ] = zero;          // all elements set to zero
  if ( Name == "" ) {
    Name = "VECTOR";
    isNameGiven = false;
  }
  isNameGiven = true;
}

static char Alphabet[] = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ";
//----------------------------------------------------------------------------------
/**
** \brief Construct a Vector from a text file
**
**    Format of the Vector file - it contains two parts: Header block and Data block.
**    HEADER BLOCK: 
**
**    (1) One or more optional COMMENT lines, where each comment line begins with '!'
**    
**    (2) Optional line with a NAME of the Vector. If given, it must be written as:
**    N= (string), 
**    where the string in the name must begin with a letter (a-w,A-W).
**    NOTE: The string must not contain spaces or tabs.
**    If the Name is not given, default name "VECTOR" is assigned.
**
**    (3) Number of elements 'n' and number of entries in a line 'l'
**    (for small vectors n and l are equal - all in one line).
**    Both entries 'n' and 'l' must be given!!
**    n= (int)  l= (int)
**           
**    The end of the Header block is given as 3 colons:
**    :::
**    DATA BLOCK:
**    numeric entris are given in a row 
**    (entries are separated by a blanko)
**
**    NOTE: For longer vectors (say, n > 15-20), multiline input is neccessary 
**    
**/
//----------------------------------------------------------------------------------
//! constructor from the file
Vector::Vector( const char* file )
{
  const int BSIZE = 120;
  FILE *from;
  int nele = 0;  // total number of elements in a vector
  int mult = 0;  // number of elements in a single line

  if ( ( from = fopen( file, "r" ) ) == 0 ) {
    cerr << "\nCannot open " << file << endl;
    exit( 1 );
  }
  char buf[ BSIZE ], *cp, *cp2;
  int efound = 0, mfound = 0, nfound = 0, colonsfound = 0;

  //--------------------------------------------- read the HEADER BLOCK
  while ( fgets( buf, BSIZE, from ) ) {
    //--- read the comment lines
    if ( ( cp = strpbrk( buf, "!" ) ) != NULL )
      * cp = '\0';
    //---- read the Name of a Vector
    if ( ( cp = strpbrk( buf, "N" ) ) != NULL )
      if ( ( cp2 = strpbrk( cp, "=" ) ) != NULL )
        if ( ( cp = strpbrk( cp2, Alphabet ) ) != NULL ) {
          Name = cp;
          isNameGiven = true;
          nfound++;
        }
    //--- read the lenth of a Vector
    if ( ( cp = strpbrk( buf, "n" ) ) != NULL )
      if ( ( cp2 = strpbrk( cp, "=" ) ) != NULL )
        if ( ( cp = strpbrk( cp2, "0123456789" ) ) != NULL ) {
          nele = atoi( cp );
          efound++;
        }
    //--- read the number of entries in a line
    if ( ( cp = strpbrk( buf, "l" ) ) != NULL )
      if ( ( cp2 = strpbrk( cp, "=" ) ) != NULL )
        if ( ( cp = strpbrk( cp2, "0123456789" ) ) != NULL ) {
          mult = atoi( cp );
          mfound++;
        }
    //--- termination of the Header block
    if ( strstr( buf, ":::" ) != NULL ) {
      colonsfound++;
      break; //... out of "while" loop
    }
  }
  if ( !efound || ! mfound || !colonsfound ) {
    cerr << "\nThe input file is not OK \n" << endl;
    exit( 1 );
  }
  if ( !nfound ) {               // no Vector name is given
    Name = "VECTOR";         // assign the default name
    isNameGiven = false;
  }
  //------------------------------------------ make space for the vector
  data = new double[ nele ];
  assert( data != NULL );
  len = nele;
  int i = 0;
  double num;
  // number of line to read
  int No_Lines, resto;
  div_t Res;
  Res.quot = 1;
  Res.rem = 0;
  if ( nele > mult )
    Res = div( nele, mult );
  No_Lines = Res.quot;        // No_Lines with full number of mult entries
  resto = Res.rem;         // ONE line with resto entries

  //-------------------------------------------- read the DATA BLOCK
  int kmax;
  for ( int j = 0; j <= No_Lines; j++ ) {   // total number: No_Lines + 1
    kmax = mult;
    if ( j == No_Lines )                     // proper number of data
      kmax = resto;                   // for the last line
    for ( int k = 0; k < kmax; k++ ) {
      char nb[ 20 ];
      fscanf( from, "%s", nb );         // scan for space-delimited string
      num = atof( nb );                 // convert string to a double value
      i = j * mult + k;               // position in the vector
      setData( i, num );                 // assign data to a vector
    }
  }
  // close the file
  fclose( from );
  if ( Name != "VECTOR" )
    isNameGiven = true;
  // clean vector of small numbers
  Clean();
}

//! copy constructor
Vector::Vector( const Vector& original )
{
  data = new double[ len = original.getLen() ];
  assert( data != NULL );
  for ( int i = 0; i < len; i++ )
    data[ i ] = original.data[ i ];
  Name = original.getName();
  if ( Name == "VECTOR" )
    isNameGiven = false;
  else
    isNameGiven = true;
}

//! copy constructor: convert from int iVector
Vector::Vector( const iVector& original )
{
  data = new double[ len = original.getLen() ];
  assert( data != NULL );
  int value;
  for ( int i = 0; i < len; i++ ) {
    value = original.getData( i );
    data[ i ] = ( double ) value;
  }
  Name = original.getName();
  if ( Name == "iVECTOR" ) {
    Name = "VECTOR";
    isNameGiven = false;
  } else
    isNameGiven = true;
}

//! destructor
Vector::~Vector()
{
  delete [] data;
}

//============================================== OPERATOR VECTOR FUNCTIONS
//
//! assignement operator
Vector& Vector::operator = ( const Vector &original )
{
  if ( this != &original ) {
    delete [] data;
    data = new double[ len = original.getLen() ];
    assert( data != NULL );
    for ( int i = 0; i < len; i++ )
      data[ i ] = original.data[ i ];
    Name = original.getName();
    if ( Name == "VECTOR" )
      isNameGiven = false;
    else
      isNameGiven = true;
  }
  return *this;
}

//=============================================== GET AND SET DATA
//
//! get an element from the Vector
const double& Vector::operator [] ( int i ) const
{
  // check index range
  assert( i>=0 && i < len );
  return data[ i ];
}

//! set an element of the Vector
double& Vector::operator [] ( int i )
{
  // check index range
  assert( i>=0 && i < len );
  return data[ i ];
}

//! get an element from the Vector
const double Vector::getData( int i ) const
{
  // check index range
  assert( i >= 0 && i < len );
  return data[ i ];
}

//! set an element to a Vector
void Vector::setData( int i, double d )
{
  // check index range
  assert( i >= 0 && i < len );
  data[ i ] = d;
}

//! set a name to a Vector
void Vector::setName( string name )
{
  Name = name;
  isNameGiven = true;
}

//===================================== ARITHMETIC OPERATORS
//
//! Vector B = self + v
Vector Vector::operator + (const Vector& v)
{
  // check vector compatibility
  assert( len == v.getLen() );
  Vector sum( len );
  for ( int i = 0; i < len; i++ )
    sum.data[ i ] = data[ i ] + v.data[ i ];
  // clean vector of small numbers
  sum.Clean();
  return ( sum );
}

//! Vector B = self - v
Vector Vector::operator - ( const Vector& v )
{
  // check vector compatibility
  assert( len == v.getLen() );
  Vector sum( len );
  for ( int i = 0; i < len; i++ )
    sum.data[ i ] = data[ i ] - v.data[ i ];
  // clean vector of small numbers
  sum.Clean();
  return ( sum );
}

//! Vector B = self * s
Vector Vector::operator * ( double s )
{
  Vector ans( len );
  for ( int i = 0; i < len; i++ )
    ans.data[ i ] = data[ i ] * s;
  // clean vector of small numbers
  ans.Clean();
  return ( ans );
}

//! scalar (dot) product of self and vector v
const double Vector::operator * ( const Vector& v ) const
{
  // check vector compatibility
  assert( len == v.getLen() );
  double s = 0.0;
  for ( int i = 0; i < len; i++ )
    s += data[ i ] * v.data[ i ];
  return s;
}

//! Vector self = self + v
Vector& Vector::operator += ( const Vector& v )
{
  // check vector compatibility
  assert( len == v.getLen() );
  for ( int i = 0; i < len; i++ )
    data[ i ] += v.data[ i ];
  // clean vector of small numbers
  Clean();
  return ( *this );
}

//! Vector self = self - v
Vector& Vector::operator -= ( const Vector& v )
{
  // check vector compatibility
  assert( len == v.getLen() );
  for ( int i = 0; i < len; i++ )
    data[ i ] -= v.data[ i ];
  // clean vector of small numbers
  Clean();
  return ( *this );
}

//! Vector self = self * s
Vector& Vector::operator *= ( double s )
{
  for ( int i = 0; i < len; i++ )
    data[ i ] *= s;
  // clean vector of small numbers
  Clean();
  return ( *this );
}

//! true if self == v
bool Vector::operator == ( const Vector& v )
{
  if ( len != v.getLen() )
    return false;
  for ( int i = 0; i < len; i++ ) {
    if ( data[ i ] != v.data[ i ] )
      return false;
  }
  if ( Name != v.getName() )    // names are not so important
    cerr << "\nVector names are NOT equal!" << endl;
  // vectors are equal
  return true;
}

//===================================== ARITHMETIC OPERATORS - equivalent functions
//
//! Vector B = self + v
Vector Vector::Add( const Vector& v )
{
  // check vector compatibility
  assert( len == v.getLen() );
  Vector sum( len );                    // mozda bi moglo: Vector sum = *this;
  for ( int i = 0; i < len; i++ )
    sum.data[ i ] = data[ i ] + v.data[ i ];  // onda bi bilo: sum.data[i] += v.data[i];
  // clean vector of small numbers
  sum.Clean();
  return ( sum );
}

//! Vector B = self - v
Vector Vector::Sub( const Vector& v )
{
  // check vector compatibility
  assert( len == v.getLen() );
  Vector sum( len );
  for ( int i = 0; i < len; i++ )
    sum.data[ i ] = data[ i ] - v.data[ i ];
  // clean vector of small numbers
  sum.Clean();
  return ( sum );
}

//! Vector B = self * s
Vector Vector::Dot( double s )
{
  Vector ans( len );
  for ( int i = 0; i < len; i++ )
    ans.data[ i ] = data[ i ] * s;
  // clean vector of small numbers
  ans.Clean();
  return ( ans );
}

//! scalar (dot) product of self and vector v
const double Vector::Dot( const Vector& v ) const
{
  // check vector compatibility
  assert( len == v.getLen() );
  double s = 0.0;
  for ( int i = 0; i < len; i++ )
    s += data[ i ] * v.data[ i ];
  return s;
}

//! outer product of self and vector v
Matrix Vector::OuterDot( const Vector& v )
{
  int m = len;
  int n = v.getLen();
  Matrix ans( m, n );
  double A_ij;
  for ( int i = 0; i < m; i++ ) {
    for ( int j = 0; j < n; j++ ) {
      A_ij = data[ i ] * v.data[ j ];
      ans.setElement( i, j, A_ij );
    }
  }
  // clean matrix of small numbers
  ans.Clean();
  return ( ans );
}

//! Add sub-vector into the vector (add to the old content)
Vector& Vector::addSubVector( int i, const Vector& Sub )
{
  assert(i >= 0);
  int r = Sub.getLen();
  // check if Sub-vector fits into this vector
  assert( (i+r) <= len );
  double ai;
  for ( int ii = 0; ii < r; ii++ ) {
    ai = this->getData( i + ii );
    this->setData( i + ii, ai + Sub.getData( ii ) );
  }
  // clean vector of small numbers
  Clean();
  return *this;
}

//! Take sub-vector from a given vector
const Vector Vector::getSubVector( int i, int subLen )
{
  assert(i >= 0 && subLen >0);
  Vector result( subLen );
  // check if subvector belongs to this vector
  assert( (i+subLen) <= len );
  double ai;
  for ( int ii = 0; ii < subLen; ii++ ) {
    ai = this->getData( i + ii );
    result.setData( ii, ai );
  }
  return result;
}

// ------------------- arithmetic operations
//! change the length of a vector
void Vector::Resize( int n )
{
  delete [] data;
  data = new double[ len = n ];
  assert( data != NULL );
}

//! set all elements to zero
void Vector::Zero( void )
{
  for ( int i = 0; i < len; i++ )
    data[ i ] = 0.0;
}

//! get max element in a vector
const double Vector::Max( void ) const
{
  double max = data[ 0 ];
  for ( int i = 1; i < len; i++ )
    if ( data[ i ] > max )
      max = data[ i ];
  return max;
}

//! get min element of a vector
const double Vector::Min( void ) const
{
  double min = data[ 0 ];
  for ( int i = 1; i < len; i++ )
    if ( data[ i ] < min )
      min = data[ i ];
  return min;
}

//! get the average value of all elements
const double Vector::Mean() const
{
  double sum = 0.0;
  for ( int i = 0; i < len; i++ )
    sum += data[ i ];
  return ( sum / len );
}

//! clean vector of small numbers 
void Vector::Clean( void )
{
  double a;
  for (int i = 0; i < len; i++) {
    a = data[ i ];
    if ( abs(a) < TINY )     // a < 1.0 e-20 
      data[ i ] = 0.0;
  }
}
 
//============================================ NORMS OF THE VECTOR
//
//! Euclid norm of the self vector
const double Vector::Norm( void ) const
{
  double s = 0.0;
  for ( int i = 0; i < len; i++ )
    s += data[ i ] * data[ i ];
  return sqrt( s );
}

//! B_Norm of the self vector: (Self_T * B * Self)^1/2
const double Vector::B_Norm( const Matrix& B ) const
{
  int n = B.getRows();
  assert( n = B.getCols() );             // matrix B must be square
  assert( len == n );
  Vector self = *this;
  Matrix B1 = B;
  Vector tmp( n );
  tmp = B1 * self;
  const double s = self.Dot(tmp);
  return sqrt( s );
}
//
//============================================= INPUT / OUTPUT FUNCTIONS
//
//----------------- console I/O
//! prints each row on a single line
void Vector::rPrint(void)
{
  cout << "Vector: " << Name << " = ";
  cout << "(";
  for ( int i = 0; i < len - 1; i++ )
    cout << data[ i ] << ",";
  cout << data[ len - 1 ] << ")" << endl;
  return ;
}

//! prints each entry on a single line
void Vector::Print(void)
{
  cout << "Vector:  " << Name << endl;
  cout << "Length = " << len << endl;
  for ( int i = 0; i < len; i++ )
    cout << data[ i ] << endl;
  return ;
}

//------------------- File I/O
//
//! read vector from predefined file
Vector Vector::Read(void)
{
  string FileName = Name + ".txt";
  Vector ans( FileName.c_str() );
  return ( ans );
}

//! read vector form a file
Vector Vector::ReadFrom( const char* file )
{
  Vector ans( file );
  return ( ans );
}

//! write vector to a predefined file
void Vector::Write(void)
{
  int L = 20;   // default number of entries in a line
  if ( len < L )
    L = len;
  string FileName = Name + ".txt";
  this->WriteTo( FileName.c_str(), L );
}

//! write vector to a file
void Vector::WriteTo( const char* file, int L )
{
  int No_Lines;          // number of lines for ONE ROW
  int resto;             // left-over in the last line for the row
  div_t Res;
  Res.quot = 1;
  Res.rem = 0;
  if ( len <= L )
    L = len;
  else if ( len > L )              // if number of cols is > number of entries in a row
    Res = div( len, L );
  No_Lines = Res.quot;        // No_Lines with full number of mult entries
  resto = Res.rem;         // ONE line with resto entries, all for just ONE ROW of data
  int J;
  string VectorName = Name;
  //-- output to file
  out.open( file, ios::out );
  assert( out != NULL );
  //--- write the HEADER BLOCK
  out << "!================== Vector(" << len << ") ====================" << endl;
  out << "N= " << VectorName << endl;
  out << "n= " << len << "  l= " << L << endl;
  out << ":::" << endl;
  //--- write the DATA BLOCK
  int kmax;
  for ( int j = 0; j <= No_Lines; j++ ) {
    kmax = L;
    if ( j == No_Lines )
      kmax = resto;
    for ( int k = 0; k < kmax; k++ ) {
      J = j * L + k;
      out << this->getData( J ) << "   ";
    }
    out << endl;
  }
  out.close();
}

//=============================================== EXTERNAL VECTOR FUNCTIONS
//
//! addition of two vectors: v3 = v1 + v2
Vector operator + ( const Vector& v1, const Vector& v2 )
{
  // check vector compatibility
  assert( v1.getLen() == v2.getLen() );
  Vector res( v1.getLen() );
  for ( int i = 0; i < v1.getLen(); i++ )
    res[ i ] = v1[ i ] + v2[ i ];
  // clean vector of small numbers
  res.Clean();
  return res;
}

//! subtraction of two vectors: v3 = v1 - v2
Vector operator - ( const Vector& v1, const Vector& v2 )
{
  // check vector compatibility
  assert( v1.getLen() == v2.getLen() );
  Vector res( v1.getLen() );
  for ( int i = 0; i < v1.getLen(); i++ )
    res[ i ] = v1[ i ] - v2[ i ];
  // clean vector of small numbers
  res.Clean();
  return res;
}

//! scalar (dot) product of two vectors:  s = v1 * v2
double operator * ( const Vector& v1, const Vector& v2 )
{
  // check vector compatibility
  assert( v1.getLen() == v2.getLen() );
  double s = 0.0;
  for ( int i = 0; i < v1.getLen(); i++ )
    s += v1[ i ] * v2[ i ];
  return s;
}

//! scalar (dot) product of two vectors: s = Dot(v1, v2)
double Dot( const Vector& v1, const Vector& v2 )
{
  // check vector compatibility
  assert( v1.getLen() == v2.getLen() );
  double s = 0.0;
  for ( int i = 0; i < v1.getLen(); i++ )
    s += v1[ i ] * v2[ i ];
  return s;
}

//! outer product of two vectors: Matrix A = v1 * v2T
Matrix OuterDot( const Vector& v1, const Vector& v2 )
{
  int m = v1.getLen();
  int n = v2.getLen();
  Matrix ans( m, n );
  double ai, bj, A_ij;
  for ( int i = 0; i < m; i++ ) {
    ai = v1.getData( i );
    for ( int j = 0; j < n; j++ ) {
      bj = v2.getData( j );
      A_ij = ai * bj;
      ans.setElement( i, j, A_ij );
    }
  }
  // clean vector of small numbers
  ans.Clean();
  return ( ans );
}

//! product of a vector and a scalar: b = s * a
Vector operator * ( double x, const Vector& v )
{
  Vector ans( v.getLen() );
  for ( int i = 0; i < v.getLen(); i++ )
    ans.setData( i, x * v[ i ] );
  // clean vector of small numbers
  ans.Clean();
  return ( ans );
}

//! product of a vector and a scalar: b = a * s
Vector operator * ( const Vector& v, double x )
{
  Vector ans( v.getLen() );
  for ( int i = 0; i < v.getLen(); i++ )
    ans.setData( i, x * v[ i ] );
  // clean vector of small numbers
  ans.Clean();
  return ( ans );
}

//! Euclid norm of a vector: s = Norm(v)
double Norm( const Vector& v )
{
  double s = 0.0;
  for ( int i = 0; i < v.getLen(); i++ )
    s += v.getData( i ) * v.getData( i );
  return sqrt( s );
}

//! B-norm (matrix norm) of a vector and a matrix: (vT * B * v)^(1/2)
double B_Norm( const Matrix& B, const Vector& v )
{
  int n = B.getRows();
  assert( n == B.getCols() );
  assert( n == v.getLen() );
  double s = v.B_Norm( B );
  return ( s );
}

//! true if two vectors are equal
bool areEqual( const Vector& v1, const Vector& v2 )
{
  //--- length of vectors
  int n = v1.getLen();
  if ( n != v2.getLen() )
    return false;
  //--- data contained in vectors
  for ( int i = 0; i < n; i++ )
    if ( v1.getData( i ) != v2.getData( i ) )
      return false;
  //--- names of Vectors are NOT important
  if ( v1.getName() != v2.getName() )
    cerr << "\nNames of vectors are NOT equal!" << endl;
  // vectors are therefore equal
  return ( true );
}

//! stream output
ostream& operator << ( ostream& s, const Vector& v )
{
  int n = v.getLen();
  s << "(";
  for ( int i = 0; i < n - 1; i++ )
    s << v.getData( i ) << ", ";
  s << v.getData( n - 1 ) << ")" << endl;
  return s;
}

//=============================================================== end of VECTOR CLASS
//
//=============================================================== MATRIX CLASS
//
//-----------------------------------------------------------------------
/**
** \brief class MATRIX of double numbers
**
**    MATRIX: MATRIX class of double numbers
**
**    Only dense matrices are considered
**    Special structure of a matrix (if any) is taken into account:
**    matrices may be diagonal, tridiagonal, symmetric or general
**/
//-----------------------------------------------------------------------
//
//! empty (default) constructor
Matrix::Matrix()
{
  m = new Vector[ M_DEFAULT ];
  assert( m != NULL );
  rows = cols = M_DEFAULT;         // Matrix (2 x 2) with all zeroes
  Vector v( cols );                  // Vector with all zeroes
  for ( int i = 0; i < rows; i++ )
    m[ i ] = v;
  Name = "MATRIX";
  Mat_Type     = NOT_CLASSIFIED;
  isSingular   = false;
  isNameGiven  = false;
  isClassified = false;
  isInverse    = false;
  isCrout      = false;
  isCholesky   = false;
}

//! Matrix with given rows, columns and all elements set to 'd'
Matrix::Matrix( int r, int c, double d )
{
  m = new Vector[ r ];
  assert( m != NULL );
  rows = r;
  cols = c;
  Vector v( cols, d );
  for ( int i = 0; i < r; i++ )  // Matrix (r x c) with all elements set to (double) d
    m[ i ] = v;
  Name = "MATRIX";
  Mat_Type     = NOT_CLASSIFIED;
  isSingular   = false;
  isNameGiven  = false;
  isClassified = false;
  isInverse    = false;
  isCrout      = false;
  isCholesky   = false;
}

//! Matrix with given rows, columns and a name (all elements set to zero)
Matrix::Matrix( int r, int c, string name )
{
  m = new Vector[ r ];
  assert( m != NULL );
  rows = r;
  cols = c;
  double d = 0.0;
  Vector v( cols, d );              // elements set to zero
  for ( int i = 0; i < r; i++ )     // Matrix (r x c) with all elements set do d
    m[ i ] = v;
  Name = name;
  if ( Name == "" ) {
    Name = "MATRIX";
    isNameGiven = false;
  }
  isNameGiven  = true;
  Mat_Type     = NOT_CLASSIFIED;
  isSingular   = false;
  isClassified = false;
  isInverse    = false;
  isCrout      = false;
  isCholesky   = false;
}
//--------------------------------------------------------- CONSTRUCT A MATRIX FROM A FILE
/**
** \brief Construct a Matrix from a file
**
**    Format of the matrix file - it contains two parts: Header block and Data block
**    HEADER BLOCK:
**
**    (1) One or more optional COMMENT lines, where each comment line begins with '!'
**    
**    (2) Optional line with a NAME of the Matrix. If given, it must be written as:
**    N= (string) 
**
**    where the string in the name must begin with a letter (a-w,A-W).
**    NOTE: The string must not contain spaces or tabs.
**    If the Name is not given, default name "MATRIX" is assigned.
**
**    (3) Number of ROWS 'r', number of COLUMNS 'c' and number of elements in one line 'l'
**    r= (int)  c= (int)  l= (int)
**
**    (for small number of columns, c = l, i.e. each row is in one line)
**    NOTE - all three entries, r, c, and l must be given!!
**
**    The end of the Header part is given as 3 colons:
**    :::
**
**    DATA BLOCK:
**    numeric entris are given row by row (number of entries in a row is <= l)
**    (entries are separated by a blanko)
**
**    NOTE: For larger matrices (say, cols > 15-20), multiline input for each row is neccessary 
**    
**/
//---------------------------------------------------------------------------------------------
//! constructor from a file
Matrix::Matrix( const char* file )
{
  const int BSIZE = 120;
  FILE *from;
  int row = 0, col = 0;       // number of rows and columns
  int mul = 0;                // number of entries in one line

  if ( ( from = fopen( file, "r" ) ) == 0 ) {
    cerr << "\nCannot open " << file << endl;
    exit( 1 );
  }
  char buf[ BSIZE ], *cp, *cp2;
  int rfound = 0, cfound = 0, mfound = 0, nfound = 0, colonsfound = 0;  // control flags

  //---------------------------------------------------- read the HEADER BLOCK
  while ( fgets( buf, BSIZE, from ) ) {
    //--- read the comment line
    if ( ( cp = strpbrk( buf, "!" ) ) != NULL )
      * cp = '\0';
    //---- read the Name of a Vector
    if ( ( cp = strpbrk( buf, "N" ) ) != NULL )
      if ( ( cp2 = strpbrk( cp, "=" ) ) != NULL )
        if ( ( cp = strpbrk( cp2, Alphabet ) ) != NULL ) {
          Name = cp;
          isNameGiven = true;
          nfound++;
        }
    //--- read the number of rows
    if ( ( cp = strpbrk( buf, "r" ) ) != NULL )
      if ( ( cp2 = strpbrk( cp, "=" ) ) != NULL )
        if ( ( cp = strpbrk( cp2, "0123456789" ) ) != NULL ) {
          row = atoi( cp );
          rfound++;
        }
    //--- read the number of columns
    if ( ( cp = strpbrk( buf, "c" ) ) != NULL )
      if ( ( cp2 = strpbrk( cp, "=" ) ) != NULL )
        if ( ( cp = strpbrk( cp2, "0123456789" ) ) != NULL ) {
          col = atoi( cp );
          cfound++;
        }
    //--- read the number of entries in a line
    if ( ( cp = strpbrk( buf, "l" ) ) != NULL )
      if ( ( cp2 = strpbrk( cp, "=" ) ) != NULL )
        if ( ( cp = strpbrk( cp2, "0123456789" ) ) != NULL ) {
          mul = atoi( cp );
          mfound++;
        }
    //--- termination of the Header block
    if ( strstr( buf, ":::" ) != NULL ) {
      colonsfound++;
      break; //... out of "while" loop
    }
  }
  if ( !rfound || !cfound || !mfound || !colonsfound ) { // if the header block is finished
    cerr << "\nThe input file is NOT OK \n" << endl;
    exit( 1 );
  }
  if ( !nfound ) {               // no Matrix name is given
    Name = "MATRIX";         // assign the default name
    isNameGiven = false;
  }
  //------------------------------------------------- make space for the matrix
  m = new Vector[ row ];             // array of row vectors
  assert ( m != NULL );
  Vector* v = new Vector( col );     // each vector with col elements
  assert ( v != NULL );
  rows = row;
  cols = col;
  int I = 0;
  double data;
  // number of input lines for each matrix row (No_Lines) and
  // number of remaining elements in the last line for a row (resto)
  // --- for one row of the matrix
  int No_Lines, resto;
  div_t Res;
  Res.quot = 1;
  Res.rem = 0;
  if ( col > mul )            // if number of cols is > number of entries in a row
    Res = div( col, mul );
  No_Lines = Res.quot;        // No_Lines with full number of mult entries
  resto = Res.rem;            // ONE line with resto entries, all for just ONE ROW of data
  //-------------------------------------- read the DATA BLOCK
  int kmax;
  for ( int i = 0; i < rows; i++ ) {
    for ( int j = 0; j <= No_Lines; j++ ) {
      kmax = mul;
      if ( j == No_Lines )                 // proper number of data
        kmax = resto;                      // for the last line in a row
      for ( int k = 0; k < kmax; k++ ) {   // read a row from a file
        char nb[ 20 ];
        fscanf( from, "%s", nb );          // scan for space-delimited string
        data = atof( nb );                 // convert string to an int value
        I = j * mul + k;
        v->setData( I, data );             // assign data to each row
      }
    }  // end of data for one row
    m[ i ] = *v;
  }
  fclose( from );
  if ( Name != "MATRIX" ) {
    isNameGiven = true;
  }
  isClassified = false;
  isSingular   = false;
  this->Classify();
  isInverse   = false;
  isCrout     = false;
  isCholesky  = false;
}

//! copy constructor
Matrix::Matrix( const Matrix& M )
{
  int i;
  rows = M.getRows();
  m = new Vector[ rows ];
  assert( m != NULL );
  cols = M.getCols();
  Vector v( cols );
  for ( i = 0; i < rows; i++ ) {
    v = M.getRowData( i );
    this->m[ i ] = v;
  }
  Name = M.getName();
  if ( Name == "MATRIX" )
    isNameGiven = false;
  else
    isNameGiven = true;
  if ( M.isClassified == true ) {
    Mat_Type = M.get_MatType();
    isClassified = true;
  } else
    this->Classify();
    isSingular = false;
    isInverse  = false;
//   isInverse = M.isInverse;
//   if (isInverse) {
//     assert ( rows == cols );
//     A_inv = new Matrix(rows,cols);
//     assert ( A_inv != NULL );
//     *A_inv = M.getInv();
//   }
  isCrout    = false;
  isCholesky = false;
}

//! copy constructor: convert from int iMatrix
Matrix::Matrix( const iMatrix& M )
{
  int i;
  rows = M.getRows();
  m = new Vector[ rows ];
  assert( m != NULL );
  cols = M.getCols();
  iVector vi( cols );
  int ni;
  double nd;
  for ( i = 0; i < rows; i++ ) {
    vi = M.getRowData( i );
    ni = vi.getData( i );
    nd = ( double ) ni;
    this->m[ i ].setData( i, nd );
  }
  Name = M.getName();
  if ( Name == "iMATRIX" ) {
    Name = "MATRIX";
    isNameGiven = false;
  } else
      isNameGiven = true;
  isSingular   = false;
  this->Classify();
  isInverse    = false;
  isCrout      = false;
  isCholesky   = false;
}

//! destructor
Matrix::~Matrix()
{
  delete [] m;
  if ( isInverse == true )
    delete A_inv;
  if ( isCrout == true )
    delete m_Crout;
  if ( isCholesky == true )
    delete m_Chol;
}

//! assignement operator
Matrix& Matrix::operator = ( const Matrix& M )
{
  if ( this != &M ) {
    delete [] m;
    rows = M.getRows();
    cols = M.getCols();
    Vector v( cols );
    m = new Vector[ rows ];
    assert ( m != NULL );
    for ( int i = 0; i < rows; i++ ) {
      v = M.getRowData( i );
      this->m[ i ] = v;
    }
    Name = M.getName();
    if ( Name == "MATRIX" )
      isNameGiven = false;
    else
      isNameGiven = true;
    Mat_Type = M.get_MatType();
    isSingular   = false;
    if ( Mat_Type != NOT_CLASSIFIED )
      isClassified = true;
    else
      this->Classify();
      isInverse = false;
//     isInverse = M.isInverse;
//     if (isInverse) {
//       assert ( rows == cols );
//       A_inv = new Matrix(rows,cols);
//       assert ( A_inv != NULL );
//       *A_inv = M.getInv();
//     }
    isCrout    = false;
    isCholesky = false;
  }
  return *this;
}
//=================================== GET ELEMENT / SET ELEMENT
//
//! get element of a matrix
const double Matrix::getElement( int i, int j ) const
{
  //check index range
  assert( i >= 0 && i < rows );
  assert( j >= 0 && j < cols );
  double val;
  val = m[ i ].getData( j );
  return val;
}

//! set a value to an element of the matrix
void Matrix::setElement( int i, int j, const double val )
{
  //check index range
  assert( i >= 0 && i < rows );
  assert( j >= 0 && j < cols );
  m[ i ].setData( j, val );
}

//! set a name to the matrix
void Matrix::setName( const string name )
{
  Name = name;
  isNameGiven = true;
}

//! get a row from a matrix
const Vector Matrix::getRowData( int i ) const
{
  //check index range
  assert( i >= 0 && i < rows );
  return m[ i ];
}

//! get a column from a matrix
const Vector Matrix::getColData( int j )const
{
  //check index range
  assert( j >= 0 && j < cols );
  int i, Rows = this->getRows();
  double val;
  Vector v(Rows);
  for ( i = 0; i < Rows; i++ ) {
    val = m[ i ].getData( j );
    v.setData( i, val );
  }
  return ( v );
}

//! set a column with a vector v
void Matrix::setColData( int j, const Vector& v )
{
  //check index range and compatibility
  assert( j >= 0 && j < cols );
  int Vrows = v.getLen();
  assert( rows == Vrows );

  int i, Rows = this->getRows();
  double val;
  for ( i = 0; i < Rows; i++ ) {
    val = v.getData( i );
    m[ i ].setData( j, val );
  }
}

//! set a row with a vector v
void Matrix::setRowData( int i, const Vector& v )
{
  //check index range and compatibility
  assert( i >= 0 && i < rows );
  int Cols = v.getLen();
  assert( cols == Cols );

  this->m[ i ] = v;
}

//! Sum of elements in a row 'i'
const double Matrix::RowSum( int i )
{
  //check index range
  assert( i >= 0 && i < rows );
  double sum = 0.0;
  for ( int k = 0; k < cols; k++ )
    sum += m[ i ][ k ];
  return sum;
}

//! Sum of elements in a column 'j'
const double Matrix::ColSum( int j )
{
  //check index range
  assert( j >= 0 && j < rows );
  double sum = 0.0;
  for ( int k = 0; k < rows; k++ )
    sum += m[ k ][ j ];
  return sum;
}

//! Matrix is assigned all elements to zero value
void Matrix::Zero( void )
{
  for ( int i = 0; i < rows; i++ )
    for ( int j = 0; j < cols; j++ )
      m[ i ][ j ] = 0.0;
}

//! Matrix is transformed into an identity matrix
void Matrix::Identity( void )
{
  // check if matrix is square
  assert ( rows == cols );
  Zero();
  for ( int i = 0; i < rows; i++ )
    m[ i ][ i ] = 1.0;
  Mat_Type = DIAGONAL;
}

//! Matrix is transformed into a diagonal matrix (all elements are equal to 'd')
void Matrix::Diagonal( double d )
{
  // check if matrix is square
  assert ( rows == cols );
  Zero();
  for ( int i = 0; i < rows; i++ )
    m[ i ][ i ] = d;
  Mat_Type = DIAGONAL;
}

//! Matrix is transformed into a diagonal matrix (elements equal to Vector 'v')
void Matrix::Diagonal( const Vector& v )
{
  // check if matrix is square
  assert ( rows == cols );
  Zero();
  for ( int i = 0; i < rows; i++ )
    m[ i ][ i ] = v.getData( i );
  Mat_Type = DIAGONAL;
}

//! find the maximum element on the main diagonal
const double Matrix::MaxDiagElem( void )
{
  // check if matrix is square
  assert ( rows == cols );
  double maxElement = getElement(0,0);
  double Aii = 0.0;
  for (int i = 1; i < rows; i++) {
    Aii = getElement(i,i);
    if (Aii > maxElement)
      maxElement = Aii;
  }
  return maxElement;
}

//! maximum element in whole matrix (square or rectangular)
const double Matrix::MaxElement( void )
{
  double maxElement, rowMax, rowEle;
  Vector Row(cols);
  maxElement = getElement(0,0);

  for (int i = 0; i < rows; i++) {
    Row = getRowData(i);
    rowMax = Row.getData(0);
    for (int j = 1; j < cols; j++) {
      rowEle = Row.getData(j);
      if (rowEle > rowMax)
        rowMax = rowEle;
   }
   if (rowMax > maxElement)
     maxElement = rowMax;
 }
 return maxElement;
}

//! clean matrix of small numbers
void Matrix::Clean( void )
{
  double Aij;
  for (int i = 0; i < rows; i++) {
    for (int j = 0; j < cols; j++) {
      Aij = getElement(i,j);
      if ( abs(Aij) < TINY )
        setElement(i,j, 0.0);
    }
  }
}

//========================================================= TRIDIAGONAL MATRICES 
//            (separate class should be developed!)
//
//! Transformation of tri-diagonal SELF matrix A(n,n) into three vectors a, b ,c
//! Vector A(n-1) contains left sub-diagonal elements of SELF matrix
//! Vector B(n)   contains diagonal elements of SELF matrix
//! Vector C(n-1) contains right sub-diagonal elements of SELF matrix
//
bool Matrix::fromTriDiag( Vector& A, Vector& B, Vector& C )
{
  // initial check
  if ( Mat_Type != TRIDIAG || Mat_Type != TRI_SYMM ) {
    cerr << "\nMatrix is NOT tri-diagonal\n";
    return false;
  }
  assert( rows == cols );  // check for square matrix
  int n = rows;
  Vector* a = new Vector( n - 1 );  // temporary working vectors
  assert( a != NULL );
  Vector* b = new Vector( n );
  assert( b != NULL );
  Vector* c = new Vector( n - 1 );
  assert( c != NULL );
  //--- select elements from tri-diagonal matrix
  int i;
  switch ( Mat_Type ) {
      case TRIDIAG:                   // non-symmetric matrix
      for ( i = 0; i < n; i++ ) {
        if ( i == 0 ) {
          b->setData( i, this->getElement( i, i ) );
          c->setData( i, this->getElement( i, i + 1 ) );
        } else if ( i < n - 1 ) {
          a->setData( i - 1, this->getElement( i, i - 1 ) );
          b->setData( i, this->getElement( i, i ) );
          c->setData( i, this->getElement( i, i + 1 ) );
        } else if ( i == n - 1 ) {
          a->setData( i - 1, this->getElement( i, i - 1 ) );
          b->setData( i, this->getElement( i, i ) );
        }
      }
      break;
      case TRI_SYMM:            // symmetric matrix: only vectors b and c, since a = c
      for ( i = 0; i < n; i++ ) {
        if ( i == 0 ) {
          b->setData( i, this->getElement( i, i ) );
          c->setData( i, this->getElement( i, i + 1 ) );
        } else if ( i < n - 1 ) {
          /*a->setData(i-1,this->getElement(i,i-1));*/
          b->setData( i, this->getElement( i, i ) );
          c->setData( i, this->getElement( i, i + 1 ) );
        } else if ( i == n - 1 ) {
          /*a->setData(i-1,this->getElement(i,i-1));*/
          b->setData( i, this->getElement( i, i ) );
        }
      }
      break;
      case NOT_CLASSIFIED:
      case DIAGONAL:
      case SYMMETRIC:
      case GEN_TYPE:
      break;
  } // end switch
  //--- now, form the final output vectors
  A.setSize( n - 1 );
  A.setName( "Left diagonal" );
  B.setSize( n );
  B.setName( "Main diagonal" );
  C.setSize( n - 1 );
  C.setName( "Right diagonal" );
  if ( Mat_Type == TRIDIAG ) {     // non-symmetric matrix
    for ( i = 0; i < n - 1; i++ ) {
      A.setData( i, a->getData( i ) );
      B.setData( i, b->getData( i ) );
      C.setData( i, c->getData( i ) );
    }
    B.setData( n - 1, b->getData( n - 1 ) );
  }
  if ( Mat_Type == TRI_SYMM ) {    // symetric matrix
    for ( i = 0; i < n - 1; i++ ) {
      A.setData( i, c->getData( i ) );   // A = C
      B.setData( i, b->getData( i ) );
      C.setData( i, c->getData( i ) );
    }
    B.setData( n - 1, b->getData( n - 1 ) );
  }
  // release memory
  delete c;
  delete b;
  delete a;

  return ( true );
}

//! Three vectors a,b,c are transformed into tri-diagonal matrix A
//! Vector b(n) is the main diagonal of tri-diagonal matrix
//! Vector a(n-1) is left and c(bn-1) is the right sub-diagonal of A
//
Matrix Matrix::setTriDiag( Vector& a, Vector& b, Vector& c )
{
  int n = b.getLen();            // order of tri-diagonal matrix
  assert( a.getLen() == n - 1 );     // check dimensions of vectors
  assert( c.getLen() == n - 1 );
  Matrix ans( n, n );
  //--- check if symmetric
  bool test = areEqual( a, c );
  if ( test == true )                 // vectors are equal: a = c
    ans.set_MatType( TRI_SYMM );
  else
    ans.set_MatType( TRIDIAG );
  //--- insert vectors into matrix
  int i;
  switch ( test ) {
      case false:           // non-symmetric tridiagonal matrix
      for ( i = 0; i < n; i++ ) {
        if ( i == 0 ) {
          ans.setElement( i, i, b.getData( i ) );
          ans.setElement( i, i + 1, c.getData( i ) );
        } else if ( i < n - 1 ) {
          ans.setElement( i, i - 1, a.getData( i - 1 ) );
          ans.setElement( i, i, b.getData( i ) );
          ans.setElement( i, i + 1, c.getData( i ) );
        } else if ( i == n - 1 ) {
          ans.setElement( i, i - 1, a.getData( i - 1 ) );
          ans.setElement( i, i, b.getData( i ) );
        }
      }
      break;
      case true:           // symmetric tridiagonal matrix
      for ( i = 0; i < n; i++ ) {
        if ( i == 0 ) {
          ans.setElement( i, i, b.getData( i ) );
          ans.setElement( i, i + 1, c.getData( i ) );
        } else if ( i < n - 1 ) {
          ans.setElement( i, i - 1, c.getData( i - 1 ) );  // a = c
          ans.setElement( i, i, b.getData( i ) );
          ans.setElement( i, i + 1, c.getData( i ) );
        } else if ( i == n - 1 ) {
          ans.setElement( i, i - 1, c.getData( i - 1 ) ); // a = c
          ans.setElement( i, i, b.getData( i ) );
        }
      }
      break;
  } // end switch test
  return ( ans );
}

//------------------------- mnipulations with subMatrices
//! Insert a sub-matrix into the matrix (replace the old content)
Matrix& Matrix::insertSubMatrix( int i, int j, const Matrix& Sub )
{
  int r = Sub.getRows();
  int c = Sub.getCols();
  // check if Sub-matrix fits into this matrix
  assert( (i+r) <= rows || (j+c) <= cols );
  for ( int ii = 0; ii < r; ii++ ) {
    for ( int jj = 0; jj < c; jj++ )
      this->setElement( i + ii, j + jj, Sub.getElement( ii, jj ) );
  }
  return *this;
}

//! Add a sub-matrix into the matrix (add to the old content)
Matrix& Matrix::addSubMatrix( int i, int j, const Matrix& Sub )
{
  int r = Sub.getRows();
  int c = Sub.getCols();
  // check if Sub-matrix fits into this matrix
  assert( (i+r) <= rows || (j+c) <= cols );
  double aij;
  for ( int ii = 0; ii < r; ii++ ) {
    for ( int jj = 0; jj < c; jj++ ) {
      aij = this->getElement( i + ii, j + jj );
      this->setElement( i + ii, j + jj, aij + Sub.getElement( ii, jj ) );
    }
  }
  // Classify();
  return *this;
}

//! Take a sub-matrix from a given matrix
Matrix Matrix::getSubMatrix( int i, int j, int subRow, int subCol )
{
  assert( subRow > 0 && subCol > 0 );
  Matrix result( subRow, subCol );
  // check if submatrix belongs to this matrix
  assert( (i+subRow) <= rows || (j+subCol) <= cols );
  double aij;
  for ( int ii = 0; ii < subRow; ii++ ) {
    for ( int jj = 0; jj < subCol; jj++ ) {
      aij = this->getElement( i + ii, j + jj );
      result.setElement( ii, jj, aij );
    }
  }
  // result.Classify();
  return result;
}

//============================================= PUBLIC FUNCTIONS: MATRIX PRODUCTS
//
//! Product of a matrix and a scalar: B = A * s (where A = SELF)
Matrix Matrix::operator * ( double s )
{
  Matrix B( rows, cols );      // result matrix: B = A * s
  B = this->Dot( s );
  return ( B );
}

//! Product of a matrix and a vector: b = A * v (where A = SELF)
Vector Matrix::operator * ( const Vector& v )
{
  // check matrix-vector comapatibility
  assert( cols == v.getLen() );
  Vector ans( rows );
  ans = this->Dot( v );
  return ( ans );
}

//! Product of matrices: C = A * B  (where A = SELF)
Matrix Matrix::operator * ( const Matrix& B )
{
  // check matrix compatibility
  assert ( cols == B.getRows() );
  int Cols = B.getCols();
  Matrix C( rows, Cols );
  C = this->Dot( B );
  return ( C );
}

//=================================== EQUIVALENT PRIVATE FUNCTIONS: MATRIX PRODUCTS
//
//! Product of a matrix and a scalar: B = A * s (where A = SELF)
Matrix Matrix::Dot( double s )
{
  Matrix B( rows, cols );      // result matrix: B = A * s
  double a_ij;
  for ( int i = 0; i < rows; i++ ) {
    for ( int j = 0; j < cols; j++ ) {
      a_ij = getElement( i, j ) * s;
      // clean matrix of small numbers
      if (abs(a_ij) < TINY)  a_ij = 0.0;
      B.setElement( i, j, a_ij );
    }
  }
  B.Classify();
  return ( B );
}

//! Product of a matrix and a vector: b = A * v (where A = SELF)
Vector Matrix::Dot( const Vector& v )
{
  // check compatibility for matrix-vector product
  assert( cols == v.getLen() );
  Vector ans( rows );               // result b
  double ai;
  for ( int i=0; i<rows; i++ ) {
    ai = 0.0;
    for ( int j=0; j<cols; j++ ) {
      ai += getElement(i,j) * v.getData(j);
      if (abs(ai) < TINY) ai = 0.0;
    }
    ans.setData( i, ai );
  }
  return ( ans );
}

//! Product of matrices: C = A * B  (where A = SELF)
Matrix Matrix::Dot( const Matrix& B )
{
  // check compatibility of matrices
  assert(cols == B.getRows());
  int Cols = B.getCols();
  Matrix C( rows, Cols );  // result matrix: C = A * B

  double Cij;
  for ( int i = 0; i < rows; i++ ) {
    for ( int j = 0; j < Cols; j++ ) {
      Cij = 0.0;
      for ( int k = 0; k < cols; k++ ) {
        Cij += getElement( i, k ) * B.getElement( k, j );
        if (abs(Cij) < TINY) Cij = 0.0;
      }
      C.setElement( i, j, Cij );
    }
  }
  C.Classify();
  return ( C );
}

//========================================= PUBLIC FUNCTIONS: ADDITION, SUBTRACTION
//
//! Addition of matrices: C = A + B  (where A = SELF)
Matrix Matrix::operator + ( const Matrix& B )
{
  // check matrix compatibility
  assert ( rows == B.getRows() && cols == B.getCols() );
  Matrix ans( rows, cols );
  ans = this->Add( B );
  return ( ans );
}

//! Substraction of matrices: C = A - B  (where A = SELF)
Matrix Matrix::operator - ( const Matrix& B )
{
  // check matrix compatibility
  assert ( rows == B.getRows() && cols == B.getCols() );
  Matrix ans( rows, cols );
  ans = this->Sub( B );
  return ( ans );
}

//! Addition of matrices: A = A + B  (where A = SELF)
Matrix& Matrix::operator += ( const Matrix& B )
{
  // check matrix compatibility
  assert ( rows == B.getRows() && cols == B.getCols() );
  *this = this->Add( B );
  return ( *this );
}

//! Substraction of matrices: A = A - B  (where A = SELF)
Matrix& Matrix::operator -= ( const Matrix& B )
{
  // check matrix comapatibility
  assert ( rows == B.getRows() && cols == B.getCols() );
  *this = this->Sub( B );
  return ( *this );
}

//================================ EQUIVALENT PRIVATE FUNCTIONS: ADDITION, SUBTRACTION
//
//! Addition of matrices: C = A + B  (where A = SELF)
Matrix Matrix::Add( const Matrix& B )
{
  // check matrix compatibility         
  // compatibility control is already done in operator functions
  assert ( rows == B.getRows() && cols == B.getCols() );
  Matrix ans( rows, cols );
  double Cij;
  for ( int i = 0; i < rows; i++ ) {
    for (int j = 0; j < cols; j++) {
      Cij = getElement(i,j) + B.getElement(i,j);
      if (abs(Cij) < TINY) Cij = 0.0;
      ans.setElement(i,j,Cij);
    }
  }
  ans.Classify();
  return ( ans );
}

//! Substraction of matrices: C = A - B  (where A = SELF)
Matrix Matrix::Sub( const Matrix& B )
{
  // check matrix compatibility
  assert ( rows == B.getRows() && cols == B.getCols() );
  Matrix ans( rows, cols );
  double Cij;
  for ( int i = 0; i < rows; i++ ) {
    for (int j = 0; j < cols; j++) {
      Cij = getElement(i,j) - B.getElement(i,j);
      if (abs(Cij) < TINY) Cij = 0.0;
      ans.setElement(i,j,Cij);
    }
  }
  ans.Classify();
  return ( ans );
}

//========================================== PUBLIC FUNCTIONS: MATRIX INVERSION
//
//! Inverse of a SELF Matrix: operator function
const Matrix Matrix::operator ! ()
{
  if ( isInverse == true )
    return ( *A_inv );
  else {
    this->Inv();
    return ( *A_inv );
  }
}

//! get inverse matrix
const Matrix Matrix::getInv()
{
  if ( isInverse == true )
    return ( *A_inv );
  else {
    this->Inv();
    return ( *A_inv );
  }
}

//================================= EQUIVALENT PRIVATE FUNCTION: MATRIX INVERSION
//
//! Inverse of a SELF Matrix
void Matrix::Inv()
{
  assert( rows == cols );                         // check if matrix is square
  A_inv = new Matrix( rows, rows );               // memory for inverse matrix
  assert( A_inv != NULL );

  //--- check if matrix is classified
  if ( Mat_Type == NOT_CLASSIFIED )
    this->Classify();

  //--- for diagonal matrix
  if ( Mat_Type == DIAGONAL && isSingular == false ) {
    double A_ii, di = 1.0;
    for ( int ii = 0; ii < rows; ii++ ) {
      A_ii = this->getElement( ii, ii );
      di *= A_ii;
      A_inv->setElement( ii, ii, 1.0 / A_ii );
    }
    Det = di;
  } // end if DIAGONAL and regular
  else if ( Mat_Type == DIAGONAL && isSingular == true ) {
    cout << "\a\nSorry, diagonal matrix is singular!" << endl;
    return;
  }  // end if DIAGONAL and singular

  //--- for symmetric matrix Cholesky decomposition
  else if ( Mat_Type == SYMMETRIC || Mat_Type == TRI_SYMM ) {
    if ( !isCholesky ) {
      m_Chol = new Cholesky( *this );
      assert( m_Chol != NULL );
      m_Chol->Decomposition();
      m_Chol->Inverse();
      *A_inv = m_Chol->getInverse();
      Det = m_Chol->getDeterminant();
      isCholesky = true;
    }
    else  if ( isCholesky ) {
      *A_inv = m_Chol->getInverse();
      Det = m_Chol->getDeterminant();
    }
  } // end if SYMMETRIC

  //--- for general matrix Crout LU decomposition
  else if ( Mat_Type == GEN_TYPE || Mat_Type == TRIDIAG ) {
    if ( !isCrout ) {
      m_Crout = new Crout( *this );
      assert( m_Crout != NULL );
      m_Crout->Decomposition();
      m_Crout->Inverse();
      *A_inv = m_Crout->getInverse();
      Det = m_Crout->getDeterminant();
      isCrout = true;
    }
    else if ( isCrout ) {
      *A_inv = m_Crout->getInverse();
      Det = m_Crout->getDeterminant();
    }
  } // end if GEN_TYPE
  string name = Name;
  name += "-Inv";
  A_inv->setName(name);
  // classify inverse matrix and set flag
  A_inv->Classify();
  isInverse = true;
  A_inv->Clean();
}

//======================================== PUBLIC FUNCTION: MATRIX TRANSPOSITION
//
//! Transpose of a SELF Matrix
const Matrix Matrix::operator ~ ()
{
  //--- check if matrix is classified
  if ( Mat_Type == NOT_CLASSIFIED )
    this->Classify();
  //--- for symmetric matrices there is no transposition
  if ( Mat_Type == DIAGONAL || Mat_Type == TRI_SYMM || Mat_Type == SYMMETRIC )
    return ( *this );
  else {                                 // for non-symmetric matrices
    Matrix ans( cols,rows );
    ans = this->Tra();
    return ( ans );         // no need to Clean(), it is already done
  }
}
//============================= EQUIVALENT PRIVATE FUNCTION: MATRIX TRANSPOSITION
//
//! Transpose of a SELF Matrix (only for non-symmetric matrices)
const Matrix Matrix::Tra( void )
{
  // matrix is not symmetric or not diagonal
  Matrix res(cols,rows);
  Vector tmp(cols);
  for (int i=0; i<rows; i++) {
    tmp = this->getRowData(i);
    res.setColData(i, tmp);
  }
  string name = Name;
  name += "-Trans";
  res.setName(name);
  res.Classify();
  return res;
}

//====================================== PUBLIC FUNCTION: MATRIX DETERMINANT
//
//! Determinant of SELF Matrix: Det = det(A)
const double Matrix::getDeterminant()
{
  if ( isInverse )              // if inverse matrix is obtained,
    return Det;                 // determinant is obtained too
  else {
    this->Inv();
    return Det;
  }
}

//====================================== PUBLIC FUNCTIONS: MATRIX DIVISION
//
//! Division of this Matrix and a Scalar: Matrix B = self / double s
Matrix Matrix::operator / ( double s )
{
  assert ( s > TINY );
  Matrix B( *this );
  double a_ij;
  for ( int i = 0; i < rows; i++ ) {
    for ( int j = 0; j < cols; j++ ) {
      a_ij = B.getElement( i, j );
      a_ij /= s;
      if ( abs(a_ij) < TINY ) a_ij = 0.0;
      B.setElement( i, j, a_ij );
    }
  }
  // clean matrix of small numbers (if s is a large number)
  // B.Clean();
  return ( B );
}

//! Division of this Matrix and a Vector
//! i.e. product A_inv * Vector v: Vector b = A_inv * v
Vector Matrix::operator / ( const Vector& v )
{
  // check compatibility
  assert ( rows == cols );                // matrix must be sqare
  assert ( rows == v.getLen() );

  Vector ans( rows );
  //--- if matrix inverse is already found
  if ( isInverse == true ) {
    ans = A_inv->Dot( v );
    // ans.Clean();  // Clean() was already done in Dot() function
    return ( ans );
  }
  else {
    this->Inv();
    ans = A_inv->Dot( v );
    // ans.Clean();
    return ( ans );
  }
}

//! Division of this Matrix and a Matrix
//! i.e. product of A_inv * Matrix B: Matrix C = A_inv * B
Matrix Matrix::operator / ( const Matrix& B )
{
  // check compatibility
  assert ( rows == cols );               // matrix must be square for inverse
  assert ( cols == B.getRows() );        // compatibility for matrix product
  // everything is OK
  int M = B.getCols();

  Matrix ans( rows, M );
  //--- if matrix inverse is already found
  if ( isInverse == true ) {
    ans = A_inv->Dot( B );
    // ans.Clean();
    return ( ans );
  }
  else {
    this->Inv();
    ans = A_inv->Dot( B );
    // ans.Clean();
    return ( ans );
  }
}

//====================================== EQUIVALENT PRIVATE FUNCTIONS: MATRIX DIVISION
//
//! Division of this Matrix and a Scalar: Matrix B = self / double s
Matrix Matrix::Div( double s )
{
  assert ( s > TINY );
  Matrix B( *this );
  double a_ij;
  for ( int i = 0; i < rows; i++ ) {
    for ( int j = 0; j < cols; j++ ) {
      a_ij = B.getElement( i, j );
      a_ij /= s;
      if ( abs(a_ij) < TINY ) a_ij = 0.0;
      B.setElement( i, j, a_ij );
    }
  }
  // clean matrix of small numbers (if s is large number)
  // B.Clean();
  return ( B );
}

//! Division of this Matrix and a Vector
//! i.e. product A_inv * Vector v: Vector b = A_inv * v
Vector Matrix::Div( const Vector& v )
{
  // check compatibility
  assert ( rows == cols );                // matrix must be sqare
  assert ( rows == v.getLen() );

  Vector ans( rows );
  //--- if matrix inverse is already found
  if ( isInverse == true )
  {
    ans = A_inv->Dot( v );
    // ans.Clean();
    return ( ans );
  }
  else {
    this->Inv();
    ans = A_inv->Dot( v );
    // ans.Clean();
    return ( ans );
  }
}

//! Division of this Matrix and a Matrix
//! i.e. product of A_inv * Matrix B: Matrix C = A_inv * B
Matrix Matrix::Div( const Matrix& B )
{
  // check compatibility
  assert ( rows == cols );               // matrix must be square for inverse
  assert ( cols == B.getRows() );        // compatibility for matrix product
  // everything is OK
  int M = B.getCols();

  Matrix ans( rows, M );
  //--- if matrix inverse is already found
  if ( isInverse == true ) {
    ans = A_inv->Dot( B );
    // ans.Clean();
    return ( ans );
  }
  else {
    this->Inv();
    ans = A_inv->Dot( B );
    // ans.Clean();
    return ( ans );
  }
}

//=========================================== PUBLIC INPUT / OUTPUT FUNCTIONS
//
//----------------- console I/O
//
//! Print each row separately (with info about the matrix)
void Matrix::Print()
{
  cout << "Matrix Name = " << Name << endl;
  cout << "No. of Rows = " << rows << " and Cols = " << cols << endl;
  cout << give_MatType() << endl;
  for ( int i = 0; i < rows; i++ )
    cout << getRowData(i);
  cout << endl;
}

//! Print each row separately (without the matrix info)
void Matrix::pPrint()
{
  for ( int i = 0; i < rows; i++ ) {
    for ( int j = 0; j < cols; j++ )
      cout << getElement(i,j) << " ";
    cout << endl;
  }
}

//------------------ File I/O
//
//!  read a Matrix from default file
Matrix Matrix::Read()
{
  string FileName = Name + ".txt";
  Matrix ans( FileName.c_str() );
  return ( ans );
}

//!  read a Matrix from a file with a given name
Matrix Matrix::ReadFrom( const char* file )
{
  Matrix ans( file );
  return ( ans );
}

//!  write Matrix to default file
void Matrix::Write()
{
  int L = 20;                 // default number of entries in a line
  if ( this->getCols() < L )
    L = this->getCols();
  string FileName = Name + ".txt";
  this->WriteTo( FileName.c_str(), L );
}

//!  write a matrix to a file with a given name
void Matrix::WriteTo( const char* file, int L )
{
  int No_Lines;          // number of lines for ONE ROW
  int resto;             // left-over in the last line for the row
  div_t Res;
  Res.quot = 1;
  Res.rem = 0;
  if ( cols <= L )
    L = cols;
  else if ( cols > L )           // if number of cols is > number of entries in a row
    Res = div( cols, L );
  No_Lines = Res.quot;      // No_Lines with full number of mult entries
  resto = Res.rem;       // ONE line with resto entries, all for just ONE ROW of data
  int J;
  double data;
  string MatrixName = Name;
  //-- output to file
  out.open( file, ios::out );
  assert( out != NULL );
  //--- write the HEADER BLOCK
  out << "!================== Matrix(" << rows << ", " << cols << ") ====================" << endl;
  out << "N= " << MatrixName << endl;
  out << "r= " << rows << "  c= " << cols << "  l= " << L << endl;
  out << ":::" << endl;
  //--- write the DATA BLOCK
  int kmax;
  for ( int i = 0; i < rows; i ++ ) {
    for ( int j = 0; j <= No_Lines; j++ ) {
      kmax = L;
      if ( j == No_Lines )
        kmax = resto;
      for ( int k = 0; k < kmax; k++ ) {
        J = j * L + k;
        data = this->getElement( i, J );
        out << data << "   ";
      }
    }
    out << endl;
  }
  out.close();
}
//=============================== PRIVATE MATRIX FUNCTIONS
//         (without equivalent public functions)
//
//! classification of a Matrix (only this is the public function)
void Matrix::Classify()
{
  //--- check if diagonal
  isDiagonal = this->isDiag();
  if ( isDiagonal == true ) {  
    Mat_Type = DIAGONAL;       
    isClassified = true;
    return ;
  }
  //--- check if tridiagonal and maybe symmetric
  int test = this->isTriDiag();
  switch ( test ) {
      case 0:
        isTridiagonal = false;
        isTriSymmetric = false;
        break;
      case 1:
        isTridiagonal = true;
        isTriSymmetric = false;
        Mat_Type = TRIDIAG;
        isClassified = true;
        break;
      case 2:
        isTriSymmetric = true;
        isTridiagonal = false;
        Mat_Type = TRI_SYMM;
        isClassified = true;
        break;
  }
  if ( isClassified == true )
    return ;
  //--- check if symmetric
  isSymmetric = this->isSymm();
  if ( isSymmetric == true ) {
    Mat_Type = SYMMETRIC;
    isClassified = true;
    return ;
  }
  //--- default case - matrix is general
  if ( isDiagonal == false && isTridiagonal == false
       && isTriSymmetric == false && isSymmetric == false ) {
    isGeneral = true;
    Mat_Type = GEN_TYPE;
    isClassified = true;
    return ;
  }
}

//! true if matrix is diagonal
bool Matrix::isDiag()
{
  // matrix must be square to be diagonal
  if ( rows != cols )
     return false;

  // check all elements of the matrix
  double Aij;
  for (int i = 0; i < rows; i++) {
    for (int j = 0; j < cols; j++) {
      Aij = getElement( i,j );
      if ( i != j && abs(Aij) > TINY )   // if at least one non-zero non-diagonal element exists,
         return false;                   // matrix is not diagonal
      if ( i == j && abs(Aij) < TINY )
         isSingular = true;              // zero element on main diagonal
    }
  }
  return true;
}

//! returns 1 if matrix is tridiagonal, 2 if tridiagonal and symmetric, 0 otherwise
int Matrix::isTriDiag()
{
  int test = 0;
  if ( rows != cols )     // matrix must be square to be tri-diagonal
    return test;
  double A_ij, A_ki;
  bool test1 = false, test2 = false;
  for ( int i = 0; i < rows; i++ ) {
    for ( int j = ( i + 1 ); j < cols; j++ ) {   
      A_ij = this->getElement( i, j );
      if ( j == ( i + 1 ) && abs( A_ij ) >= TINY )
        test1 = true;
      if ( j > ( i + 1 ) && abs( A_ij ) >= TINY )
        return test;
    }
    for ( int k = ( i + 1 ); k < cols; k++ ) {
      A_ki = this->getElement( k, i );
      if ( k == ( i + 1 ) && abs( A_ki ) >= TINY )
        test2 = true;
      if ( k > ( i + 1 ) && abs( A_ki ) >= TINY )
        return test;
    }
  }
  if ( test1 == true && test2 == true ) {     // matrix is tri-diagonal
    // now, check if symmetric
    for ( int i = 1; i < rows; i++ ) {
      A_ij = this->getElement( i - 1, i );    // right sub-diagonal
      A_ki = this->getElement( i, i - 1 );    // left sub-diagonal
      if ( A_ij != A_ki )                     // if only one pair is .ne.
        return ( test = 1 );                  // matrix is not symmetric
    }
    test = 2;                                 // matrix is symmetric and tridiagonal
  }
  return test;
}

//! true if matrix is symmetric
bool Matrix::isSymm()
{
  if ( rows != cols )                         // matrix must be square to be symmetric
    return false;
  double A_ij, A_ji;
  for ( int i = 0; i < rows - 1; i++ ) {
    for ( int j = ( i + 1 ); j < cols; j++ ) {
      A_ij = this->getElement( i, j );
      A_ji = this->getElement( j, i );
      if ( A_ij != A_ji )                   // if at least one pair of symmetric elements are not equal
        return false;                       // matrix is NOT symmetric
    }
  }
  return true;  // matrix is symmetric
}

//! give the matrix type as a String
string Matrix::give_MatType()  const
{
  char * base = "Matrix is ";
  string info( base );
  switch ( Mat_Type )
  {
      case NOT_CLASSIFIED:
      info = info + "NOT_CLASSIFIED\n";
      break;
      case DIAGONAL:
      info = info + "DIAGONAL\n";
      break;
      case TRIDIAG:
      info = info + "TRIDIAGONAL\n";
      break;
      case TRI_SYMM:
      info = info + "TRIDIAGONAL SYMMETRIC\n";
      break;
      case SYMMETRIC:
      info = info + "SYMMETRIC\n";
      break;
      case GEN_TYPE:
      info = info + "GENERAL\n";
      break;
  }
  return ( info );
}

//=============================== EXTERNAL MATRIX FUNCTIONS
//
//!  Multiply a Matrix with a Scalar
Matrix Dot( double x, Matrix& M )
{
  Matrix ans( M.getRows(), M.getCols() );
  double Aij, Mij;
  for ( int i = 0; i < ans.getRows(); i++ ) {
    for ( int j = 0; j < ans.getCols(); j++ ) {
      Mij = M.getElement( i, j );
      Aij = Mij * x;
      if ( abs(Aij) < TINY ) Aij = 0.0;
      ans.setElement( i, j, Aij );
    }
  }
  // clean matrix of small numbers
  // ans.Clean();
  return ( ans );
}

//!  Multiply a Matrix with a Scalar
Matrix Dot( Matrix& M, double x )
{
  Matrix ans( M.getRows(), M.getCols() );
  double Aij, Mij;
  for ( int i = 0; i < ans.getRows(); i++ ) {
    for ( int j = 0; j < ans.getCols(); j++ ) {
      Mij = M.getElement( i, j );
      Aij = Mij * x;
      if ( abs(Aij) < TINY ) Aij = 0.0;
      ans.setElement( i, j, Aij );
    }
  }
  // clean matrix of small numbers
  // ans.Clean();
  return ( ans );
}

//! Multiplication of two matrices: C = A * B
Matrix Dot( Matrix& A, Matrix& B )
{
  assert( A.getCols() == B.getRows() );
  int M = A.getRows();
  int N = B.getCols();
  Matrix C(M,N);
  C = A * B;
  // clean matrix of small numbers
  // C.Clean();
  return C;
}

//============================================== LU decomposition for tri-diagonal matrices
/**
** \brief LU decomposition for tri-diagonal matrices: A = L * U
**
**  Tri-diagonal matrices may be symmetric or non-symmetric
**  INPUT:  tri-diagonal matrix A(n,n) - destroyed after lu_triDiag
**  OUTPUT: LU decomposition -  A = L * U
**  
**  Matrix U has 'ones' as diagonal elements and upper diagonal
**  Matrix L has diagonal elements and left diagonal
**  Therefore, matrices L*U are presented as one tri-diagonam matrix
**  contained in input matrix A
**
**  After LU decomposition, determinent, inverse matrix or 
**  solution of algebraic equations may be obtained
**/
//---------------------------------------------------------------------------------------
void lu_triDiag( Matrix& A )
{
  int n = A.getRows();
  assert( n == A.getCols() );              // check if A is square matrix

  Vector* a = new Vector( n - 1 );         // temporary working vectors
  assert( a != NULL );                     // to represent the matrix A
  Vector* b = new Vector( n );
  assert( b != NULL );
  Vector* c = new Vector( n - 1 );
  assert( c != NULL );
  //--- get diagonal vectors
  bool test = A.fromTriDiag( *a, *b, *c );
  if ( test == false ) {
    cerr << "\nProblems with transformation to three vectors\n ";
    return ;
  }

  Vector* Alfa = new Vector( n - 1 );  // temporary working vectors
  assert( Alfa != NULL );              // to make the LU tridiagonal matrix
  Vector* Beta = new Vector( n );
  assert( Beta != NULL );
  Vector* Gama = new Vector( n - 1 );
  assert( Gama != NULL );

  int i;
  double Ai, Bi, Ci;
  double Di, Ei;
  for ( i = 0; i < n; i++ ) {
    if ( i == 0 ) {
      Bi = b->getData( i );
      Ci = c->getData( i );
      Beta->setData( i, Bi );
      Gama->setData( i, Ci / Bi );
    } else if ( i < n - 1 ) {
      Ai = a->getData( i );
      Bi = b->getData( i );
      Ci = c->getData( i );
      Alfa->setData( i, Ai );
      Di = Ai * Gama->getData( i - 1 );
      Beta->setData( i, Bi - Di );
      Ei = Ci / ( Bi - Di );
      Gama->setData( i, Ei );
    } else if ( i == n - 1 ) {
      Ai = a->getData( i );
      Bi = b->getData( i );
      Alfa->setData( i, Ai );
      Di = Ai * Gama->getData( i - 1 );
      Beta->setData( i, Bi - Di );
    }
  } // end for i
  //--- now, construct tridiagonal LU matrix
  A.setTriDiag( *Alfa, *Beta, *Gama );

  // release dynamic memory
  delete Gama;
  delete Beta;
  delete Alfa;
  delete c;
  delete b;
  delete a;
}

//============================ Inverse of a square Matrix A (remains the same)
//
//! Inverse of a square Matrix A
Matrix Inv( Matrix& A )
{
  Matrix * tmp = new Matrix( A );       // copy of the input Matrix A
  assert( tmp != NULL );
  int N = A.getRows();
  assert ( N == A.getCols() );         // check if Matrix A is square

  Matrix ans( N, N );
  ans = ! ( *tmp );
  // release dynamic memory
  delete tmp;
  // ans.Clean();  // already done Clean() in Inv
  // give solution
  return ( ans );
}

//============================== Transpose of a Matrix A (remains the same)
//
//! Transpose of a Matrix A
Matrix Tra( Matrix& A )
{
  Matrix * tmp = new Matrix( A );       // copy of the input Matrix A
  assert( tmp != NULL );
  int N = A.getRows();
  int M = A.getCols();                 // SELF may be N x M MAtrix

  Matrix ans( M, N );
  ans = ~( *tmp );
  // release dynamic memory
  delete tmp;
  // give solution
  return ( ans );
}

//=============================== Determinant of a Matrix A
//
//! Determinant of a Matrix A
double Det( Matrix& A )
{
  int N = A.getRows();
  assert ( N == A.getCols() );
  if ( A.isInverse == true )
    return A.getDeterminant();
    
  Matrix* tmp = new Matrix( A );
  assert( tmp != NULL );
  tmp->getInv();
  double D = tmp->getDeterminant();
  delete tmp;
  return D;
}

//! Stream output
//
ostream& operator << ( ostream& s, const Matrix& M )
{
  s << M.give_MatType();
  for ( int i = 0; i < M.getRows(); i++ )
    s << M.getRowData(i);
  return s;
}

}   // end of namespace smv_lib
//============================================================== end of MATRIX CLASS
//============================================= end of file
